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u Abstract 

In this paper we report on the feasibility of light confinement in orbital geodesies on station- 
ary, planar, and centro- symmetric refractive index mappings. Constrained to fabrication and 

i [meta] material limitations, the refractive index, n, has been bounded to the range: 0.8 < n (r) < 3.5. 
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trappings by generalized orbit conditions defined a priori. Our simulation results show that the 



above mentioned refractive index distributions trap light in an open orbit manifold, both perennial 
and temporal, in regards to initial conditions. Moreover, due to their characteristics, these map- 
pings could be advantageous to optical computing and telecommunications, for example, providing 



an on-demand time delay or optical memories. Furthermore, beyond their practical applications to 
photonics, these mappings set forth an attractive realm to construct a panoply of celestial mechanics 



analogies and experiments in the laboratory. 
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I. INTRODUCTION 

Einstein's general relativity theory allows us to understand the dynamics of mass and 
massless particles as they travel through space-time. From it, we learn that matter and 
geometry are inextricably related; the latter transfigures the metric of space-time, and this 
alteration in turn modifies the geodesies, the paths of particles, resulting in a manifold 
of relativistic effects, v. gr. the apparent shift in the position of the stars due to the 
gravitational lensing effect, depicted in Fig. [Taj where photons traveling near a massive 
star are deflected from their otherwise straight path due to the curvature of space-time; or 
trapping in multifarious orbits around black holes with respect to the initial conditions and 



shape of the gravitational attractor [IH3], see Figs, lb and lc 

Keeping in mind the previous illustrations, let us momentarily examine a more mundane 
experiment. Imagine a photon moving in a region of space filled with some monotonically 
increasing refractive index. What we observe as the photon moves in the medium is that its 
path describes a curve like the one depicted in Fig. [TdJ If the gradient of the refractive index 



is tailored properly, it would be possible to build an optical analogy to gravitational bending 



(see Fig. la) for example. This leads us to wonder whether it is possible to fabricate an 
optical analogue to other sidereal geodesies by means of a differential control on the index 
of refraction. 

The latter analogy reveals a connection between both examples, a weft where both realms 
[sidereal and laboratory] can be understood as manifestations of the same mathematical 
ground, albeit significantly different physics, differential geometry. Einstein's general theory 
of relativity tells us that geometry and gravity are linked by means of the metric tensor 
of space-time; in optics the intricacy of the relation between geometry and the physical 
properties of materials (permittivity and permeability) has been evinced by miscellaneous 
methods, most based on differential geometry [IHTj. 

These methods can be used to show that the optical example of differential light bending 
can be understood in terms of geometry, essentially by bending the space on which the light 
propagates, hence the metric is related to the electromagnetic tensors [HE]. This realization 
paved the way for the theoretical realization of optical cloaks, super-lenses, optical illusion 
devices and recently black-hole analogues in moving refractive index perturbations [6l-fT4]. 

Concerning trapping and light confinement in particular, different methods have been 




(a) 



(b) 




S^6 



:^ggfe 







(c) 



7 ; >9^^- 



(d) 



FIG. 1: Comparison of sidereal and laboratory frame geodesies for light beams in (a) the 

vicinity of a star, (b) and (c) spinning black hole at different initial conditions, and (d) the 

laboratory across an heterogeneous monotonic increasing refractive index medium. 



proposed [2l [T5UT7] . Of particular interest is Genov et al study, which presented an isotropic 
inhomogeneous refractive index mapping which was able to confine light mimicking a closed 
orbit around a black hole [15J. As the authors noted, such mapping could be implemented 
on a dielectric photonic crystal, composed of air-gaps in Gallium- Arsenide. 

In our work, we posit theoretical traps that provide an optical analogue to celestial 
mechanics to realize light confinement, trapping in a region of space, by means of two- 
dimensional stationary refractive index mappings, which could be implemented under cur- 
rent technological and photo-refractive [met a] materials, constraint to optical frequencies, 
which could be easily integrable on all-optical-chips [T8H2T] . The mathematical and physi- 
cal background to make these effects possible could bring forth an exciting ground to test 
sidereal mechanics in the laboratory; and provide the key to enable a miscellany planar 
optical systems that are of great interest for diverse photonic applications, ranging from 
time delays, to temporal optical memories, or random resonators, to name a few; thence its 
significance. 



II. DERIVING THE MODEL 

We regard a model equation based in the paraxial approximation of non-linear electrody- 
namics, by looking into a medium where the index of refraction can be modulated in space 
and time by some smooth varying function /, e.g. in photo-refractive media index of refrac- 
tion is a real valued function / of the incident light intensity /; therefore n{x) = /(/(x)), 
where x is the space-time quadrivector which corresponds to the coordinates {x° : x 1 , x 2 , x 3 }. 
In this representation x° corresponds to the time coordinate, and the rest to the usual three- 
dimensional space. In this context the photon's path is governed by the geometry of the 
space the light is propagating on, which in general has the metric tensor <jy, and the in- 
finitesimal length element ds 2 = gijdx l dxK Therefore, if we parameterized the path followed 
by light with reference to the proper time r, i.e. x(t), we can analytically derive the path 
equations by means of Lagrange-Euler equations, whit a Lagrangian of the form: 

dx l <9x J 
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(1) 



Or dr 
where the derivatives are taken in respect to the affine parameter r, i.e. the proper time. 

The physical description behind Eq. ([!]) is relevant to optics, as we set forth earlier in the 
introduction. If the metric were indeed related to the electromagnetic tensors, permittivity 
and permeability, then it could be possible to write the light-path in the paraxial approx- 
imation as a function of either the refractive index n = /(x°,r) = yjs{r^ t)/x(r, £), or the 
metric ds. Hence, it could be possible to use Fermat's principle, where Eq. ([!]) would take 
the role of an optical Lagrangian, i.e. the solution to the variational problem which states 
that given a refractive index n(x°, r) the path light follows is given by the extrema of: 

Sr 



= sfds = 0, (2) 



X\ 



where the integration is done between the arbitrary points x\ to x 3 2) with j = 0, . . . , 3, and 5 
indicates the functional maximization. Recursively, the solution to Eq. ^ is the Lagrange- 
Euler Eq. ([l]). Meanwhile, observe that in the integral of Eq. ^ the term ds takes into 
account the metric by virtue of the covariant metric tensor g^ of the curved space: 

ds = ^/gijdx l dxi . (3) 

Yet, we are still to determine the relation between the metric tensor g^ and the refractive 
index n. This can be achieved by means of conformal mappings, where it is possible to build 



the relation between the metric and the electromagnetic properties of a [meta] material [IH6]. 
Specifically, if the line element of space-time takes the form ds 2 = g odx dx — gudx l dx\ it 
is possible to map the curved space-time to an isotropic, non-dispersive and non-absorbing 
medium as [H [15] : 

—g.gv 

n = e = fi= a , (4) 

V #00 

where g = det (gij). Furthermore if we consider Minkowski's metric of special relativity, the 
tensor g has the entries: 

g = diag(g 00l g lu #22, #33) = diag(l, -1, -1, -1), (5) 

then the refractive index n = n(x°,r) can be found using Eq. Q. And consequently, the 
length element ds of Eq. ^ can be written as: 

ds 2 = —dx°dx° - dx l dx\ i = 1, 2, 3 (6) 

A. Path equations and the inverse problem 

Accordingly, from the Lagrange-Euler (fTl) we can derive the light-path, which in cylin- 
drical coordinates takes the form: 

t d Tt 

r = r(/) 2 -{ *z— , (7) 
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(9) 
t ( td t n + 2 ( zd z n + <pdj>n + rd r n ) ) 

i=-± ^ y~, (10) 

n 
where all derivatives marked by □ and are with respect to the affine parameter r, e.g. 
f = dr/dr. Equation ^ through Eq. ^ describe the parametrized path of a light beam 
traveling in a space-time with such geometry or equivalent refractive index. 

We will describe light being trapped, when the geodesies are perpetually bounded to 
some region of space, for as long as the material properties remain unchanged. In this sense 
a light trapping geodesic is one where any of the following conditions arise: 

1. r,0 -» 0, 



2. r — >• and m ^ < 6 < 



mm — Y — Ymaxi 



3. -► and r min < r < r max , 



4. r min < r < r marr and 6 min <6< 



for all r > r*; where r* marks the moment at which the initial conditions for trapping are 
met. The first condition is only possible if the refractive index n goes to infinity, or if the 
refractive index co-moves with the wave [HI H7j . The other three cases can be realized in 
stationary refractive index [T5l [22]. Thus, they can be imposed as general conditions to 
solve the inverse problem to Eq. ^ through Eq. (10). 



The inverse problem begins by defining a set of a priori conditions that describes a 
generalized trapping scenario. If we impose premises [T] to [4j described earlier, the path 
equations reduce to new set of equations where the affine parameter is time itself; whence 
we can rework Eq. ([7|) through Eq. (10) into polar coordinates: 



nrS 2 — 2r0cLn — f 2 d r n + r 2 6 2 d r n , ^ N 

r= ^ , 11 
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Since we are interested in the general stationary case, we impose the fourth condition 



above, ensuing < r < r max and (j) m in < < </Wr, on Eq. (11) and Eq. (12), and find a 
non-singular refractive index distribution, n(r, 0), that satisfies them. Notice that we require 
the refractive index to be smoothly varying and C 2 continuous, i.e. Vn • e 1 < M, where M is 
a real valued number, and e 1 a vector in a give direction in space-time. The general solution 
to this problem, however, is complex and a unique solution is not guaranteed |T8l E3J l2lj . 

The latter being said, the problem can be simplified further if we recall that a trapping 
orbit occurs when the radius r is bounded to a region of space, i.e. r min < r < r maX) for 
all 0; whence r(r max ) = i(r min ) = 0, and r can be written in terms of the angle (/) or the 
time t. Hence a natural solution to the inverse problem is an autonomous refractive index 
that varies with the radial variable, i.e. n = n(r). This refractive index map resembles the 
centro-symmetric characteristics of gravitational fields, where, as pointed out independently 
by Genov and Ni pL5ll22] . light would behave similarly to that traveling in curved space-time. 
Such distribution can be constrained to the conditions listed previously. 

This result is significant in terms of fabrication where novel photo-refractive [metal] materials 
could be used. The refractive index in these composites can be modulated by the intensity 
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of the change-inducing beam; which can be controlled by shaping its transversal profile 
[T8H2T] . Accordingly several families of refractive index distributions can be achieved, e.g. 
Gaussian, Laguerre, or Bessel, which are all centro-symmetrical. 



B. System stability and dynamic equillibrium 



Since the refractive index is a function of the radius, Eq. (11) and Eq. (12) can be 



rewritten in the form (13): 



x = —h{x,x) — g{x) ) 

r — ( r 2 (j) 2 — r 2 ) d r In n + rcj) 2 , 

•• . • n + rd r n 

= — zrcp , 

nr 



(13) 
(14) 

(15) 



the system's symmetry now allows us to integrate Eq. (15), hence: 
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where £ is a constant of integration. This equation is reminiscent of orbital mechanics, 
with the sole exception that /x, the mass in celestial mechanics, varies with the radius. 



Substituting the latter result into Eq. (14), and rearranging yields: 

02 
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(17) 



Equation (17) arrangement allows us to calculate the ratio of external energy supply, 

S]; simplifying we get: 



which is a measure of the equilibrium stability of the system 



St = — = —fh(r, f)—f 



n 4 r 2 



r d r lnn. 



(18) 



By virtue of Bertrand's theorem [T5l [27], we know that f ^0, for all time t; then from 
premises [3] and [4] it follows that r min < r < r maX) where r min < 0; and from the definition 
of n(r) as a monotonous decreasing function it follows that d r lnn < 0. Therefore, as the 
radius decreases, the radial velocity goes to its minimum and St < 0, therefore the system 
approaches the equilibrium point and the amplitude of the path decreases. On the other 
hand, as the radius increases, the velocity tends to a maximum, and S t > 0, thus the system 
is driven away from the equilibrium point and the amplitude of the path increases. In Fig. 
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FIG. 2: i?aizo o/ external energy supply, 8 t) for a radial symmetric system, where x = r 

and y — f. Notice that there exists a x c , critical radius, at wich for all phase space pairs 

{x, y\x > x c } £ t = 0, which implies that the light geodesies do not describe an orbital 

motion. 



[2] we show £ t for a Gaussian like refractive index described in section |IH A[ Observe that in 
the light confinement scenario we will see an oscillatory behavior in £ t . 

To conclude we study the stability of the orbits, by means of Lyapunov's stability 



theory for the governing Eq. (17), and derive the constraints applicable to the refractive 



index mapping to attain stable orbital paths [15]. First, let us rewrite Eq. (17) as: 



r = x, 



k= (r 2 f3 2 -H 2 ) X (r)+r(3\ 



(19) 



with (3 = <j> and x(r) = d r lnn. Then, if we introduce small perturbations r = r + 5 and 
k = x + e, expand %(r) in a Taylor series around r , and drop all higher order terms for S 



and e; we can rewrite Eq. (19) as a linear system of equations for the perturbation terms: 



V 




(20) 



x oX o, W = 2 X o(3 2 r + x'o(r 2 o(3-<), 



where Xo an d x'o are the h near terms of the Taylor expansion of %(r) around r . Ensuing 
the eigenvalues of the matrix above are given by: 



X 12 = V± W 2 + W. 



(21) 



The geodesies will be stable if the real part of the eigenvalues of Eq. ( |20| ), ^{Xi^}, are 
negative or equal to zero. Recall that < f3 < lim r ^ 00 n(r) _1 , and x(r) < 0; and therefore 
the condition for Lyapunov stability is reached if: 

d 



dr X(r) 



x < 0, V + W < 0, (22) 



This constrained is met for all points r near the stability point r = 0, provided the 
radial velocity meets r < 0. In general, since the refractive index is monotonous decrescent, 
the system will be Lyapunov stable at points where r < 0, see, for example, the phase-space 
plots in Figs. [6b] and [8b} In addition, it is interesting to analyze whether the system will 
be Lyapunov stable under the conditions necessary to sustain a circular orbit. In this case 
k = 0, and r = a, constant, for all time t] consequently V = 0, W < 0, and the second term 



in the equality (22) is imaginary, provided that x( r ) ^ 0; thus resulting in ^{Xi^} = V = 0. 
These conditions, however, result in a different family of refractive indexes. In our study 
we focus on Gaussian-like refractive index distributions, feasible through photo-refraction, 
where x/0, ergo the system will not describe closed circular orbits. 

In summary, light confinement will occur as long as the conditions cited earlier apply. 
Under them the system will oscillate around the equilibrium point as depicted by the ratio 
of external energy supply. In the next section we will present the confinement geodesies 
resulting from these refractive index distributions. We will show that the geodesies in this 
system mimic those of a spinning black hole and other celestial objects, thus reproducing 
the paths of light traveling in curved space. 

III. RESULTS AND DISCUSSION 

In this section we explore three refractive index distributions based on Gaussian profiles. 
The first two mimic light trapped in the vicinity of a black hole, and the third case that of 
a particle moving in-between two planets. 

A. Guassian distribution, single attractor 

A Gaussian refractive index distribution has the form: 

~ r2/a2 -u - (23) 
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FIG. 3: Light path dynamics for Gaussian refractive index shape for different n a and n c . 

The maximum value of n in all cases is 3.8 at the center of the device, a is constant for all 

results. Notice the increment in the internal radius r min as n c increases, and the orbit 

width reduces. 

Then making it possible to constrain the maximum and minimum of the refractive index 
distribution to feasible values, e.g. 0.8n < 3.8. Hence, given a range of values for the 
background and maximum refractive index, n c and n a , respectively, we can evaluate different 
light confinement scenarios, see Fig. [3| 

As discussed in the previous section the orbital paths are open, while their shape and 
stability depend significantly on the initial conditions, n(r) and {r , r , O , O }. As portrayed 
by Fig. [3] different initial conditions will result in new orbit manifolds with distinct charac- 
teristics, albeit maintaining a common morphology. As it can be seen from these figures the 
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orbital paths are highly sensitive to the modifications in the refractive index and the initial 
conditions; a change in ~ 0.5% in n c alters the orbit, as seen from the plot of radii extrema 
in Fig. |U 
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FIG. 4: Inner and outer radius as a function of the background refractive index n c . As it 
increases we observe that the outer radius diminishes, while the inner increases. 
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(a) (b) 

FIG. 5: Confinement orbit for light traveling in a Gaussian attractor where n a = 3.0 and 

n c = 0.8. (a) Shows a schematic representation of the planar trapping device where the 

refractive index increases towards the center, where red indicates maximum refractive index 

n = 3.8, and light blue n = 0.8; (b) polar-2D path as calculated by the simulation routine. 

In Fig. [5] we show the 2D open orbit for a refractive index map where n c = 0.8 and 



n a = 3.0. As detailed in Fig. [6aJ the radial velocity describes an oscillatory form, which 
indicates that the radius follows a similar motion, swinging between two extrema, as expected 
in a confine and open orbital motion. The latter can be observed in greater detail by 
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FIG. 6: Phase space for Gaussian-like refractive index map with n a = 3.0, n c = 0.8. The 
plots show: (a) radial r and angular cp velocities, (b) radial phase space, and (c) radial 

velocity as a function of the angle (c). 



analyzing the phase-space, plotted in Fig. 6b; notice that when the velocity r reaches zero 
the orbit radius r reaches an extrema; for a closed single loop the amplitude of the velocity 
oscillation would be zero, and the phase space would become a horizontal line at r — 0. The 
behaviour of the angle versus the radial velocity, and consequently to the radius movement, 



is described by Fig. [6cJ where the absolute value of the radial velocity is plotted against 
the angular position, which results in a periodic motion. Observe that in the limit at every 
angle there would be a time at which the radial velocity would be zero, which is typical of 
open orbits. 



B. Mexican hat distribution 

We study a variation on the Gaussian trap to understand the potential modifications on 
the trapping orbits when large modifications ~ 20% of the refractive index take place. The 
mathematical description of the mexican hat is: 

r2 



n d r _ r 2 /a 2 
n= {n a — )e ' 



+ n c 



(24) 



Using the same conditions as in the previous analysis we set n c = 0.8, n a = 3.0, and set 
rid = 0.2. The results, shown in Figs. [7]and|8j display an orbit whose appearance resembles 
that of the regular Gaussian trap studied earlier. However, a closer look reveals that the trap 
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confinement area is <~ 2% slimmer, and the oscillatory frequency of the radial and angular 



velocities, Fig. [8aJ increases by ~ 30%. Whereas the radial phase-space, Fig. [8bJ shows 
a smaller eccentricity compared to its Gaussian counterpart, and the absolute value of the 



radial velocity versus the angle, Fig. 8c, describes the typical behaviour of an open orbit. 
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FIG. 7: (a) Schematic of Mexican hat like refractive index map where n a = 3.0 and 

n c = 0.8 and depression is n^ = 0.2; (b) depicts the resulting trapped orbit in polar 

coordinates as calculated by the simulation routine. 

These variations in the trapping conditions could potentially be used to differentiate 
trapping in transient optical memories to allocated binary data. 

C. Double attractors 

Interested in observing phenomenologically trapping in a binary system we evinced a 
planar refractive index map resulting from the superposition of two Gaussian functions, 



n = n a (e-( r - r °ff l)2/(j2 + e -(*-w/2) 2 /* 2 ) + ^ 



(25) 



The resulting orbits are shown in Fig. [9| A binary system, however, presents a relatively 
higher sensitivity to modifications in the refractive index, where a distortion of An > 0.03 
modifies the trapping orbits, in some cases releasing the light beam. Yet, it is possible to 
confine the light beam to an orbit which mimics the behaviour of Newton's three body 
problem, where, as expected, small perturbations result in sundry possible solutions to the 
path equation, i.e. chaos. 
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FIG. 8: Phase-space for a mexican hat refractive index map withn a = 3.0, n c = 0.8. The 
plots show: (a) radial r and angular (j) velocities, (b) radial phase space, and (c) radial 

velocity as a function of the angle (c). 




(a) n a = 2.7, n c = 1.1 



(b) n a = 2.9, n c = .9 



FIG. 9: Geodesies for light in a binary system. Observe that trapping also occurs for 
certain initial conditions at specific background and maximum refractive index. 

D. The space-time equations, trapping 



Briefly, we show that the same Gaussian solution explored in section |HI A| fulfills the 
trapping conditions in the full eikonal space-time, Eq. ([7]) through Eq. ([9]). The results are 



shown in Fig. 10, 

In this case we observe a similar behaviour to that of the Gaussian trap in the stationary 
approximation, the difference in the amplitude values is due to normalization. An oscillatory 
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FIG. 10: Phase space of the trapping orbit for n a = 3.0, n c = 0.8. The plots show: (a) 
radial r and angular cp velocities, (b) radial phase space, and (c) radial velocity as a 

function of the angle. 

behaviour in the radius and angular velocities is described in Figs. 10(b) and 10(c); as it 
would be expected from a centro-symmetrical refractive index mapping. 

E. Material realization 



Recent advancement in photo-refractive materials [T8H2T] could enable the construction 
of the light confinement devices discussed. For example, a prototype Gaussian trap could 
be fabricated by stacking different photo-refractive [meta] materials like As x Si- x in multiple 
layers. Recalling that in photo-refraction the optical properties of the material, permeability 
and permittivity, change as a function of the light intensity, allowing the fabrication of a 
continuously varying refractive index as the one required by a centro-symmetrical attractors. 
Moreover, since photo-refractivity is elastic, i.e. reversible, it allows one to modify the 
refractive index of the layers combined in time, thence enabling or disabling the trap at will. 



IV. CONCLUSIONS 

We have demonstrated that, under the paraxial description of light in heterogeneous 
media, it is possible to use centro-symmetrical refractive index distributions as attractors 
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for light, which are able to perform light trapping or confinement in open orbits. These 
mappings are bound within fabrication constrains, and due to the paraxial approximation 
they are several times larger than the mean beam width and wavelength, ^ 50 times larger. 
The proposed mappings have potential applications to transient optical memories, delays, 
concentrators, random cavities, and beam stirrers, to cite a few; which could enable the next 
generation of photonic integrated circuits, PICs, processors and computation systems. Due 
to their sensitivity to the change in the refractive index these devices could be employed as a 
sensor platform. Furthermore, since they are drawn in close analogy to celestial phenomena 
they could be used as laboratory setups to test a large variety of effects in celestial mechanics. 
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